function [Y,X] = simu_1st(gx, hx, eta, sig, e)

%xp^i = hx^i_a x_a  + sig * eta^i_c ep_c
%
%y^i = gx^i_a x_a 
%
%where 
% hx is nx by nx
% gx is ny by nx
% eta is nx by ne
% sig is a parameter scaling the innovation to the exogenous variables 
% e is T by ne random shock
% T is the length of the simulation 
%
%Output:  vectors X (T by nx) and Y (T by ny) containing time series for x and y and e (T by 1) containing  
%  modified
% (c) Stephanie Schmitt-Grohe and Martin Uribe, January 22, 2002

T=size(e,1);
nx=size(hx,1);
ny=size(gx,1);

X=zeros(T,nx);
Y=zeros(T,ny);

X(1,:)=(sig*eta*e(1,:)')';
Y(1,:)=(gx*X(1,:)')';
for i=2:T;
 X(i,:)=(hx*X(i-1,:)'+sig*eta*e(i,:)')';
 Y(i,:)=(gx*X(i,:)')';
end





